OHI British Columbia | OHI Science | Citation policy
OHI British Columbia | OHI Science | Citation policy
knitr::opts_chunk$set(fig.width = 6, fig.height = 4, fig.path = 'figs/',
echo = TRUE, message = FALSE, warning = FALSE)
# library(ggmap)
library(rgdal)
library(sf)
library(broom)
source('~/github/ohibc/src/R/common.R') ### an OHIBC specific version of common.R
dir_git <- '~/github/ohibc'
dir_spatial <- file.path(dir_git, 'prep/_spatial') ### github: general buffer region shapefiles
dir_anx <- file.path(dir_M, 'git-annex/bcprep')
### goal specific folders and info
goal <- 'ao'
scenario <- 'v2017'
dir_goal <- file.path(dir_git, 'prep', goal, scenario)
dir_goal_anx <- file.path(dir_anx, goal, scenario)
### provenance tracking
library(provRmd); prov_setup()
### set up the default BC projection to be BC Albers
p4s_bcalb <- c('bcalb' = '+init=epsg:3005')This data is incorporated into the OHI British Columbia First Nations Resource Access and Opportunities (AO) goal.
Fishing licenses
Fishery Management Area boundaries
Census subdivision administrative boundaries
2016 Census Total Population Results Census Subdivisions
Read in DFO contamination closures from .xlsx files.
lic_clean_file <- file.path(dir_goal_anx, 'int/licenses_clean.csv')
license_file <- file.path(dir_anx, '_raw_data/fishing_license/Pacific Region Commercial Fishing Licenses En.csv')
metadata_file <- file.path(dir_anx, '_raw_data/fishing_license/Commercial License Meta Data EN (1).csv')
if(!file.exists(lic_clean_file)) {
licenses_raw <- read_csv(license_file) %>%
setNames(tolower(names(.)) %>% str_replace_all('[\\s()]+', '_'))
for(tmp in names(licenses_raw)) {
### tmpname <- names(x)[1]
licenses_raw[[tmp]] <- ifelse(is.na(licenses_raw[[tmp]]) & class(licenses_raw[[tmp]]) == 'character',
'',
licenses_raw[[tmp]])
}
metadata <- read_csv(metadata_file, skip = 54) %>%
setNames(tolower(names(.)) %>% str_replace_all(' ', '_'))
licenses_clean <- licenses_raw %>%
left_join(metadata, by = c('licence_prefix' = 'prefix')) %>%
mutate(lic_holder = ifelse(licence_holder_type == 'VESSEL',
paste('vessel', licence_holder_vessel),
paste(licence_holder_first_name, licence_holder_midddle_name, licence_holder_last_name)),
lic_oper = ifelse(licence_holder_type == 'VESSEL',
paste('vessel', licence_operator_vessel),
paste(licence_operator_first_name, licence_operator_midddle_name, licence_operator_last_name)),
vessel_own = ifelse(licence_holder_type == 'VESSEL',
paste(vessel_owner_first_name, vessel_owner_midddle_name, vessel_owner_last_name),
NA)) %>%
mutate(year = licence_suffix,
fn = str_detect(tolower(prefix_description), 'aborig|northern native|northernnative'),
fn = ifelse(is.na(fn), str_detect(tolower(licence_holder_type), 'aborig'), fn)) %>%
select(year,
lic_prefix = licence_prefix,
prefix_desc = prefix_description,
lic_tab = licence_tab,
fn, area,
fishery_group, foreign_vessel,
lic_hold_type = licence_holder_type,
lic_hold_vessel = licence_holder_vessel,
lic_oper_vessel = licence_operator_vessel,
lic_holder, lic_oper, vessel_own,
vessel_length_m = vessel_length_mr_)
write_csv(licenses_clean, lic_clean_file)
} else {
message('File exists: ', lic_clean_file)
git_prov(c(license_file, metadata_file), 'input')
git_prov(lic_clean_file, 'output')
}This link provides a look-up of different commercial management areas to PFMAs and PFMSAs. This might be useful for spatializing the licenses.
http://www.pac.dfo-mpo.gc.ca/fm-gp/licence-permis/areas-secteurs-eng.html
The table in this link has been saved as a .csv in the ao/v2017/raw folder.
This bar chart shows the proportion of licenses allocated to First Nations for each fishery group. The indicated proportion is the maximum for that fishery group across all years.
lic_clean <- read_csv(lic_clean_file)
license_fn_yr <- lic_clean %>%
select(year,
prefix_desc,
fn,
lic_prefix,
fishery_group) %>%
mutate(fishery_group = tolower(fishery_group)) %>%
group_by(year, fishery_group) %>%
summarize(n_total = n(),
n_fn = sum(fn),
pct_fn = n_fn / n_total) %>%
ungroup()
license_fn_yr <- license_fn_yr %>%
group_by(fishery_group) %>%
mutate(max_pct_fn = max(pct_fn, na.rm = TRUE)) %>%
ungroup()
bar_by_gp <- license_fn_yr %>%
select(max_pct_fn, fishery_group) %>%
filter(max_pct_fn > 0) %>%
distinct() %>%
arrange(desc(max_pct_fn)) %>%
mutate(fishery_group = factor(fishery_group, levels = unique(fishery_group))) %>%
ggplot(aes(x = fishery_group, y = max_pct_fn)) +
geom_bar(stat = 'identity') +
theme(axis.text.x = element_text(angle = 75, hjust = 1))
print(bar_by_gp)
This examines the total fishing license allocation to First Nations across all license types, fishery groups, and regions, by year.
lic_clean <- read_csv(lic_clean_file)
license_fn_yr1 <- lic_clean %>%
select(year,
prefix_desc,
fn,
lic_prefix,
fishery_group) %>%
mutate(fishery_group = tolower(fishery_group)) %>%
group_by(year) %>%
summarize(n_total = n(),
n_fn = sum(fn),
pct_fn = n_fn / n_total) %>%
ungroup()
lic_total_plot <- ggplot(license_fn_yr1, aes(x = year, y = n_total)) +
geom_area(fill = 'grey60', color = 'grey40', size = .25) +
geom_area(aes(y = n_fn), fill = 'steelblue3', color = 'grey40', size = .25) +
labs(y = 'Number of licenses')
lic_pct_plot <- ggplot(license_fn_yr1, aes(x = year, y = pct_fn)) +
geom_line(color = 'steelblue3', size = .5) +
labs(y = 'Percent FN licenses')
print(lic_total_plot)
print(lic_pct_plot)
Much of this code was developed for the AO Closures prep script; since there are likely idiosyncracies within the data here compared to the closures data, the code has been copied and customized rather than functionalized.
Clean up area and subarea calls from the cleaned data. No need to keep the license info at this point; can rejoin later. Many of the descriptions are complex text such as “Area 12, Area 13, except SA 13-1 to 13-12, 13-15 to 13-17, Area 14 (except SA 14-5, 14-8, 14-15, Area 15 (except SA 15-5), Areas 16 to 29 (except SA 17-20, 23-6 and 29-5)” which must be parsed down to areas.
lic_area_all <- read_csv(lic_clean_file) %>%
select(area, year) %>%
group_by(area) %>%
mutate(year_1 = min(year, na.rm = TRUE),
year_n = max(year, na.rm = TRUE)) %>%
ungroup() %>%
select(-year) %>%
distinct() %>%
mutate(area = ifelse(is.na(area), 'no area designated', tolower(area)))
lic_area_lookup <- read_csv(file.path(dir_goal, 'raw/license_area_lookup.csv')) %>%
mutate(lic_area = tolower(lic_area),
desc = tolower(desc))
# x <- lic_area_all %>%
# filter(area %in% lic_area_lookup$lic_area) %>%
# arrange(area)
# y <- lic_area_all %>%
# filter(!area %in% lic_area_lookup$lic_area) %>%
# arrange(area)
# z <- lic_area_lookup %>%
# filter(!lic_area %in% lic_area_all$area) %>%
# arrange(lic_area)
### All the remaining unmatched license areas (y) seem to have ended long ago,
### most before 2000, and only a couple as late as 2002 or 2003.
### Presumably these license areas are no longer valid, so don't show up
### in the license area designations.
### Because we are counting within regions rather than area-weighting, it is
### less important to worry about subareas. Drop them for simplicity. This
### also means not having to worry about "except" or "excludes"
lic_lookup_clean <- lic_area_lookup %>%
select(lic_area, desc) %>%
mutate(area_edit = desc,
area_edit = str_replace_all(area_edit, '(?<=[0-9])[-—.](?=[0-9])', '~'), ### this line writes subareas as ##~##
area_edit = str_replace_all(area_edit, '\\s+', ' '))
### step 2: expand 'to' for areas (e.g. areas 3 to 10).
### * Identify clauses, split into 'from'/'to'
### * create a vector of from:to and convert to string.
### * str_replace the clause with the expanded vector.
### * reassemble the entire closure list string
lic_lookup_to <- lic_lookup_clean %>%
filter(str_detect(area_edit, '[0-9] to [0-9]')) %>%
mutate(area_split = str_split(tolower(area_edit), ',|;|and')) %>% ### divide at comma or semicolon
unnest(area_split) %>%
mutate(to_clause = str_extract(area_split, '[0-9]*~?[0-9]+ to [0-9]+~?[0-9]*')) %>%
separate(to_clause, into = c('from_num', 'to_num'), sep = '[^0-9]? to ', remove = FALSE) %>%
rowwise() %>%
### build vector of Area numerals, skipping sub-areas
mutate(to_vec = ifelse(!is.na(from_num) & !str_detect(from_num, '~') & !str_detect(to_num, '~'),
paste(as.integer(from_num):as.integer(to_num), collapse = ', '),
'')) %>%
### build vector of Sub-Area numerals, skipping areas
mutate(tmp_area = ifelse(!is.na(from_num) & str_detect(from_num, '~') & str_detect(to_num, '~'),
str_extract(from_num, '[0-9]+(?=~)'), NA),
to_vec_sub = ifelse(!is.na(from_num) & str_detect(from_num, '~') & str_detect(to_num, '~'),
paste(tmp_area, as.integer(str_extract(from_num, '(?<=~)[0-9]+')):as.integer(str_extract(to_num, '(?<=~)[0-9]+')), sep = '~', collapse = ', '),
'')) %>%
ungroup() %>%
mutate(to_vec = ifelse(to_vec == '', to_vec_sub, to_vec),
area_split2 = ifelse(!is.na(to_clause),
str_replace(area_split, to_clause, to_vec),
area_split),
area_split3 = str_replace_all(area_split2, '[^0-9~]+', ' '),
area_split3 = str_trim(area_split3)) %>%
select(-area_split, -area_split2, -to_clause, -from_num, -to_num, -to_vec) %>%
distinct() %>%
filter(!str_detect(area_split3, '[0-9]{4}')) %>% ### ditches four digit #s (i.e. years)
group_by(lic_area) %>%
summarize(desc = first(desc),
area_edit = paste(area_split3, collapse = ' ')) %>%
ungroup()
lic_lookup_all <- lic_lookup_clean %>%
filter(!lic_area %in% lic_lookup_to$lic_area) %>%
mutate(area_edit = str_replace_all(area_edit, '[0-9]{4}', ''),
area_edit = str_replace_all(area_edit, '[^0-9~]+', ' '),
area_edit = str_trim(area_edit)) %>%
bind_rows(lic_lookup_to) %>%
mutate(pfma_area = str_split(area_edit, ' ')) %>%
unnest(pfma_area) %>%
separate(pfma_area, c('pfma', 'pfmsa'), '~') %>%
mutate(pfma = as.integer(pfma),
pfmsa = as.integer(pfmsa)) %>%
select(-area_edit) %>%
distinct()
### attach OHIBC regions by PFMA
pfma_lookup <- read_csv(file.path(dir_git, 'prep/fis/v2017/int/pfmsa_to_ohibc.csv')) %>%
select(rgn_id, pfma_id, pfmsa_id) %>%
filter(!is.na(pfma_id) & !is.na(rgn_id)) %>%
distinct()
lic_lookup_rgn <- lic_lookup_all %>%
filter(is.na(pfmsa)) %>%
inner_join(pfma_lookup, by = c('pfma' = 'pfma_id')) %>%
bind_rows(lic_lookup_all %>%
filter(!is.na(pfmsa)) %>%
inner_join(pfma_lookup, by = c('pfma' = 'pfma_id', 'pfmsa' = 'pfmsa_id'))) %>%
select(-pfmsa_id, -pfmsa) %>%
distinct()
write_csv(lic_lookup_rgn, file.path(dir_goal, 'int/lic_lookup_rgn.csv'))lic_lookup_rgn <- read_csv(file.path(dir_goal, 'int/lic_lookup_rgn.csv')) %>%
filter(!is.na(rgn_id))
lic_clean <- read_csv(lic_clean_file) %>%
mutate(area = tolower(area))
lic_fn_by_rgn <- lic_clean %>%
inner_join(lic_lookup_rgn, by = c('area' = 'lic_area')) %>%
select(-pfma) %>%
distinct()
x <- lic_clean %>%
group_by(area) %>%
mutate(year_1 = min(year, na.rm = TRUE),
year_n = max(year, na.rm = TRUE)) %>%
ungroup()
y <- x %>%
filter(!year <= 1995) %>%
filter(!str_detect(area, 'yukon|stikine|taku')) %>%
filter(!is.na(area))
### still some areas not bound; these all end pre-2003, most pre-1999
sum_fn_by_rgn <- lic_fn_by_rgn %>%
group_by(year, rgn_id) %>%
summarize(n_fn = sum(fn),
n_tot = n(),
pct_fn = n_fn/n_tot) %>%
filter(year >= 1995) %>% ### no First Nations licenses before 1995 on record
filter(year < 2017) ### incomplete current year
write_csv(sum_fn_by_rgn, file.path(dir_goal, 'output/ao_licenses.csv'))sum_fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
left_join(get_rgn_names(), by = 'rgn_id') %>%
mutate(pct_fn = 100 * pct_fn)
plot_fn_by_rgn <- ggplot(sum_fn_by_rgn, aes(x = year, y = pct_fn, group = rgn_name, color = rgn_name)) +
ggtheme_plot() +
geom_line(size = 1.5, alpha = 0.7) +
labs(y = 'Percent of licenses to First Nations',
color = 'OHIBC Region') +
scale_color_brewer(palette = 'Dark2')
print(plot_fn_by_rgn)
bc_poly_84 <- sf::read_sf(dsn = dir_spatial %>% path.expand(),
layer = 'ohibc_rgn_simple') %>%
st_transform(4326) %>%
mutate(rgn_id = as.integer(rgn_id))
sum_fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
mutate(pct_fn = pct_fn * 100) %>%
filter(year %% 5 == 0)
bc_poly_fn <- bc_poly_84 %>%
full_join(sum_fn_by_rgn, by = 'rgn_id')
ao_max_pct <- bc_poly_fn$pct_fn %>% max()
bpMap <- ggplot(bc_poly_fn) +
ggtheme_plot() +
theme(axis.text = element_blank(),
axis.title = element_blank()) +
geom_sf(aes(fill = pct_fn),
color = 'slateblue', size = .25, alpha = .4) +
scale_fill_gradientn(colours = brewer.pal(7, 'RdYlBu'), limits = c(0, ao_max_pct)) +
labs(title="OHIBC AO",
fill = '% licenses\nFirst Nation') +
facet_wrap( ~ year)
ggsave(filename = file.path(dir_goal, 'int/ao_licenses_mapped.png'),
plot = bpMap)
One possible reference point is that the % of licenses apportioned to FN should be proportional to the % of population in region that is First Nations.
ohibc_rgn <- read_sf(dsn = dir_spatial,
layer = 'ohibc_rgns_unclipped')
csd_clean <- function(csd_sf, year) {
csd_sf2 <- csd_sf %>%
clean_df_names() %>%
select(csduid, csdname, csdtype, geometry) %>%
st_transform(st_crs(ohibc_rgn))
csd_sf2$a_tot <- st_area(csd_sf)
csd_sf_bc <- csd_sf2 %>%
st_buffer(dist = 0) %>%
st_intersection(ohibc_rgn)
csd_sf_bc$a_ohibc <- st_area(csd_sf_bc)
csd_bc_df <- csd_sf_bc %>%
as.data.frame() %>%
clean_df_names() %>%
mutate(prop_area = as.numeric(a_ohibc / a_tot),
area_total_km2 = as.numeric(a_tot / 1e6),
area_ohibc_km2 = as.numeric(a_ohibc / 1e6),
year = year) %>%
select(csduid, csdname, csdtype, rgn_id,
area_total_km2, area_ohibc_km2, prop_area, year)
return(csd_bc_df)
}
if(!file.exists(file.path(dir_goal, 'int/ao_licenses_csd_2016.csv'))) {
csd_2016 <- read_sf(dsn = file.path(dir_anx, '_raw_data/bcstats/csd2016carto'),
layer = 'lcsd000b16a_e') %>%
csd_clean(2016)
write_csv(csd_2016, file.path(dir_goal, 'int/ao_licenses_csd_2016.csv'))
}pop_2016 <- read_csv(file.path(dir_anx, '_raw_data/bcstats',
'2016 Census - CSD_control_79d4a796-2ccd-4732-a01c-1a8fd26d254a_1.csv')) %>%
clean_df_names() %>%
select(sgc, name, type,
pop = `2016_population`,
prev_pop = `2011_population_2`)
write_csv(pop_2016, file.path(dir_goal, 'int/ao_licenses_pop_2016.csv'))pop_df <- read_csv(file.path(dir_goal, 'int/ao_licenses_pop_2016.csv'))
csd_df <- read_csv(file.path(dir_goal, 'int/ao_licenses_csd_2016.csv'))
pop_csd <- pop_df %>%
left_join(csd_df, by = c('sgc' = 'csduid')) %>%
filter(!is.na(rgn_id)) %>%
select(-type, -name, -prev_pop) %>%
group_by(rgn_id, csdtype) %>%
summarize(pop = sum(pop * prop_area, na.rm = TRUE))
write_csv(pop_csd, file.path(dir_goal, 'int/ao_licenses_pop_csd_2016.csv'))Develop a dataframe of municipalities by OHIBC region (inland, based on MaPP inland regions and watershed boundaries). Census subdivisions will be allocated to OHIBC regions according to area-weighting.
First Nations communities can be identified from the census data using the “type” codes:
| CSD type | Description |
|---|---|
| CY | City |
| DM | District Municipality |
| T | Town |
| VL | Village |
| IM | Island Municipality |
| IRI ✓ | Indian Reserve |
| RDA | Regional District Electoral Area |
| IGD ✓ | Indian Governmental District |
| S-E ✓ | Indian Settlement |
| NL ✓ | Nisga’a Land |
| NVL ✓ | Nisga’a Village |
Note, this does not catch First Nations members who are living separately from these communities, and counts all members of the communities as First Nations members.
pop_csd <- read_csv(file.path(dir_goal, 'int/ao_licenses_pop_csd_2016.csv'))
fn_types <- c('IRI', 'IGD', 'S-E', 'NL', 'NVL') %>%
paste(collapse = '|')
pop_fn_rgn <- pop_csd %>%
group_by(rgn_id) %>%
mutate(fn_pop = str_detect(csdtype, fn_types)) %>%
summarize(fn_pop = sum(pop * fn_pop),
tot_pop = sum(pop),
pct_fn_pop = round(fn_pop / tot_pop, 5))
knitr::kable(pop_fn_rgn, caption = 'First Nations communities as percent of total pop')| rgn_id | fn_pop | tot_pop | pct_fn_pop |
|---|---|---|---|
| 1 | 3178.346 | 4.415680e+04 | 0.07198 |
| 2 | 1435.355 | 4.328024e+03 | 0.33164 |
| 3 | 2259.224 | 3.625213e+03 | 0.62320 |
| 4 | 2336.893 | 3.324516e+04 | 0.07029 |
| 5 | 4937.957 | 3.328312e+05 | 0.01484 |
| 6 | 28014.948 | 3.253829e+06 | 0.00861 |
| 8 | 0.000 | 1.009863e+00 | 0.00000 |
write_csv(pop_fn_rgn, file.path(dir_goal, 'output/ao_licenses_fn_pop.csv'))license_ref_min <- 0.20 ### arbitrary lower limit of reference point
license_status <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
filter(rgn_id != 7) %>%
left_join(read_csv(file.path(dir_goal, 'output/ao_licenses_fn_pop.csv')),
by = 'rgn_id') %>%
mutate(ref_pt = ifelse(pct_fn_pop > license_ref_min, pct_fn_pop, license_ref_min),
status = pct_fn / ref_pt,
status = ifelse(status > 1, 1, status),
component = 'first_nations_licenses') %>%
select(year, rgn_id, pct_fn, pct_fn_pop, ref_pt, status, component) %>%
left_join(get_rgn_names(), by = 'rgn_id')
plot_fn_by_rgn <- ggplot(license_status, aes(x = year, y = pct_fn, group = rgn_name, color = rgn_name)) +
ggtheme_plot() +
geom_hline(aes(yintercept = pct_fn_pop), color = 'grey', size = 1.5, alpha = .5) +
geom_hline(aes(yintercept = ref_pt), color = 'red', alpha = 1) +
geom_line(size = 1.5, alpha = 0.7) +
labs(y = 'Percent of licenses to First Nations',
color = 'OHIBC Region') +
scale_color_brewer(palette = 'Dark2') +
facet_wrap( ~ rgn_name)
print(plot_fn_by_rgn)
Note: grey line indicates % of rgn population that are members of First Nations communities according to 2016 census subdivisions; red line indicates reference point as proportion of region population that is First Nations, with a minimum set at 20%.
pop_fn_rgn <- read_csv(file.path(dir_goal, 'int/ao_licenses_fn_pop.csv'))
fn_by_rgn <- read_csv(file.path(dir_goal, 'output/ao_licenses.csv')) %>%
left_join(get_rgn_names(), by = 'rgn_id') %>%
left_join(pop_fn_rgn, by = 'rgn_id') %>%
filter(rgn_id != 7) %>%
rename(pct_fn_lic = pct_fn)
fn_by_rgn_mean <- fn_by_rgn %>%
group_by(rgn_id, rgn_name, pct_fn_pop, fn_pop) %>%
summarize(pct_fn_lic = mean(pct_fn_lic, na.rm = TRUE),
n_fn = mean(n_fn, na.rm = TRUE)) %>%
ungroup()
x <- ggplot(fn_by_rgn, aes(x = pct_fn_pop, y = pct_fn_lic)) +
geom_abline(slope = 1, intercept = 0, color = 'red') +
geom_point(aes(color = rgn_name), alpha = .3) +
geom_point(data = fn_by_rgn_mean, aes(color = rgn_name), size = 4) +
labs(title = 'FN % licenses vs % rgn pop')
y_mdl <- lm(n_fn ~ fn_pop, data = fn_by_rgn)
y_coef <- c('slope' = y_mdl$coefficients[['fn_pop']],
'intcp' = y_mdl$coefficients[['(Intercept)']],
'r2' = y_mdl %>% summary() %>% .$r.squared)
mean_ratio <- fn_by_rgn %>%
group_by(rgn_id) %>%
summarize(mean_n = mean(n_fn),
mean_fn_pop = mean(fn_pop),
mean_ratio = mean(n_fn / fn_pop)) %>%
ungroup() %>%
filter(rgn_id != 8) %>%
mutate(unif_dist = sum(mean_n) / sum(mean_fn_pop))
y <- ggplot(fn_by_rgn, aes(x = fn_pop, y = n_fn)) +
geom_smooth(method = 'lm', color = 'red') +
geom_point(aes(color = rgn_name), alpha = .3) +
geom_point(data = fn_by_rgn_mean, aes(color = rgn_name), size = 4) +
labs(title = 'FN licenses vs rgn pop')
print(x); print(y)Do licensed vessels differ greatly between the First Nations fleet and the non-First Nations fleet? A significant difference in boat size could indicate a significant difference in the access to catch for one license type vs. another.
boats <- read_csv(lic_clean_file, nogit = TRUE) %>%
select(year, fishery_group, fn, length = vessel_length_m, lic_holder, lic_oper) %>%
mutate(length = ifelse(length == 0, NA, length))
bigboats <- boats %>%
filter(length > 50)
medboats <- boats %>%
group_by(year, fn) %>%
mutate(mean_l = mean(length, na.rm = TRUE)) %>%
ungroup() %>%
filter(length <= 40) %>%
filter(year %in% c(1995, 2000, 2005, 2010, 2015))
length_hist <- ggplot(medboats, aes(x = length, fill = fn)) +
geom_histogram(aes(frame = year), alpha = .6) +
geom_vline(aes(xintercept = mean_l, color = fn)) +
geom_text(data = medboats %>%
select(year, mean_l, fn) %>%
distinct(),
aes(x = mean_l, y = 1, label = paste0('mean = ', round(mean_l, 1), ' m')),
hjust = 0, vjust = 0, angle = 90, size = 3) +
labs(x = 'vessel length, meters',
y = 'license count',
fill = 'First Nations?',
color = 'First Nations?') +
facet_wrap( ~ year)
# print(length_hist)
ggsave(file.path(dir_goal, 'int/ao_licenses_vessel_lengths.png'))
lic_lookup_rgn <- read_csv(file.path(dir_goal, 'int/lic_lookup_rgn.csv')) %>%
filter(!is.na(rgn_id))
lic_clean <- read_csv(file.path(dir_goal_anx, 'int/licenses_clean.csv')) %>%
mutate(area = tolower(area))
lic_by_pfma_all <- lic_clean %>%
inner_join(lic_lookup_rgn, by = c('area' = 'lic_area')) %>%
select(-rgn_id) %>%
distinct()
### Total licenses by fishery group:
# lic_by_pfma_all$fishery_group %>% table()
# CLAM CRAB GEODUCK AND HORSE CLAM HERRING HERRING - FOOD & BAIT HERRING - SPECIAL ISSUE
# 113384 48600 4024 74 4760 348
# OYSTERS, PACIFIC PROCESSING RED SEA URCHIN ROCKFISH ROE HERRING GILL NET ROE HERRING SEINE
# 1650 283 14776 331645 89560 16515
# SALMON GILL NET SALMON HISTORICAL SALMON SEINE SALMON TROLL SEA CUCUMBER
# 316457 20 97852 208194 6047
### FN licenses by fishery group:
# lic_by_pfma_all %>% filter(fn) %>% .$fishery_group %>% table()
# CLAM CRAB GEODUCK AND HORSE CLAM OYSTERS, PACIFIC RED SEA URCHIN ROCKFISH ROE HERRING GILL NET
# 24405 1898 9 156 1266 12061 3750
# ROE HERRING SEINE SALMON GILL NET SALMON SEINE SALMON TROLL SEA CUCUMBER
# 236 75222 4271 6914 10
### These have NO fn licenses:
# lic_by_pfma_all %>% group_by(fishery_group) %>% filter(sum(fn) == 0) %>% .$fishery_group %>% table()
# HERRING HERRING - FOOD & BAIT HERRING - SPECIAL ISSUE PROCESSING SALMON HISTORICAL
# 74 4760 348 283 20
lic_by_pfma <- lic_by_pfma_all %>%
select(year, lic_tab, fn, fishery_group, desc, pfma) %>%
mutate(taxa = case_when(str_detect(fishery_group, 'CLAM|OYSTER') ~ 'bivalve',
str_detect(fishery_group, 'SALMON') ~ 'salmon',
str_detect(fishery_group, 'HERRING') ~ 'herring',
str_detect(fishery_group, 'CRAB') ~ 'crab',
str_detect(fishery_group, 'URCH|CUCUM') ~ 'echinoderm',
str_detect(fishery_group, 'ROCKFISH') ~ 'rockfish',
str_detect(fishery_group, 'PROCESS') ~ 'processing',
TRUE ~ 'other')) %>%
group_by(year, pfma, taxa, fn) %>%
summarize(n_lic = n()) %>%
ungroup()
write_csv(lic_by_pfma, file.path(dir_goal, 'int/lic_by_pfma.csv'))
lic_by_year <- lic_by_pfma %>%
group_by(year, taxa) %>%
summarize(fn_lic = sum(n_lic * fn),
tot_lic = sum(n_lic)) %>%
ungroup() %>%
filter(year >= 2000 & year <= 2016) %>%
filter(taxa != 'processing')
plot_by_year <- ggplot(lic_by_year, aes(x = year)) +
ggtheme_plot() +
geom_line(aes(y = fn_lic, group = taxa),
color = 'darkred', alpha = .7, size = 1) +
geom_line(aes(y = tot_lic, group = taxa),
color = 'darkblue', alpha = .5, size = 1) +
scale_color_brewer(palette = 'Dark2') +
facet_wrap( ~ taxa, scales = 'free_y') +
ylim(c(0, NA))
print(plot_by_year)
ggsave(file.path(dir_goal, 'figs/licenses/lic_by_year_and_taxa.png'))lic_by_pfma_tot <- read_csv(file.path(dir_goal, 'int/lic_by_pfma.csv')) %>%
filter(year >= 2000 & year <= 2016) %>%
group_by(pfma, taxa) %>%
summarize(fn_lic = sum(n_lic * fn),
tot_lic = sum(n_lic)) %>%
ungroup()
pfma_sf <- sf::read_sf(dsn = file.path(dir_anx, '_raw_data/dfo_khunter',
'management_boundaries/d2016',
'pac_fishery_mgmt_areas'),
layer = 'DFO_BC_PFMA_AREAS_50K_V3_1') %>%
select(pfma = STATAREA)
bbx <- st_bbox(pfma_sf)[c(1, 3, 2, 4)] %>%
raster::extent()
garbage <- capture.output({
continent_sf <- readOGR(dsn = path.expand(file.path(dir_spatial)),
layer = 'ohibc_continent') %>%
raster::crop(bbx) %>%
sf::st_as_sf()
})
pfma_lic_sf <- pfma_sf %>%
left_join(lic_by_pfma_tot, by = 'pfma') %>%
filter(!is.na(taxa))
for(taxa_sel in unique(pfma_lic_sf$taxa)) {
### taxa_sel <- 'herring'
# cat(taxa_sel, '\n')
tmp_df <- pfma_lic_sf %>%
filter(taxa == taxa_sel) %>%
mutate(tot_lic = ifelse(tot_lic == 0, NA, tot_lic),
fn_lic = ifelse(fn_lic == 0, NA, fn_lic))
breaks <- c(1, 10, 100, 1000, 10000)
max_tot <- tmp_df$tot_lic %>% max(na.rm = TRUE)
breaks_tot <- breaks[breaks < max_tot * 10]
tot_map <- ggplot(tmp_df %>% filter(!is.na(tot_lic))) +
ggtheme_plot() +
geom_sf(data = continent_sf, fill = '#bb9999', color = NA) +
geom_sf(data = pfma_sf, fill = 'grey96',
color = 'black', size = .1) +
geom_sf(aes(fill = tot_lic),
color = NA, alpha = .8) +
scale_fill_distiller(palette = 'Blues',
breaks = breaks_tot,
limits = c(1, NA),
trans = 'log', direction = 1) +
labs(title = taxa_sel,
fill = 'Total licenses')
print(tot_map)
ggsave(file.path(dir_goal, sprintf('figs/licenses/lic_map_pfma_%s_total.png', taxa_sel)))
tmp_df_fn <- tmp_df %>% filter(!is.na(fn_lic))
if(nrow(tmp_df_fn) > 0) {
max_fn <- tmp_df$fn_lic %>% max(na.rm = TRUE)
breaks_fn <- breaks[breaks < max_fn * 10]
fn_map <- ggplot(tmp_df_fn) +
ggtheme_plot() +
geom_sf(data = continent_sf, fill = '#bb9999', color = NA) +
geom_sf(data = pfma_sf, fill = 'grey96',
color = 'black', size = .1) +
geom_sf(aes(fill = fn_lic),
color = NA, alpha = .8) +
scale_fill_distiller(palette = 'Blues',
breaks = breaks_fn,
limits = c(1, NA),
trans = 'log', direction = 1) +
labs(title = taxa_sel,
fill = 'FN licenses')
print(fn_map)
ggsave(file.path(dir_goal, sprintf('figs/licenses/lic_map_pfma_%s_fn.png', taxa_sel)))
}
}
prov_wrapup(commit_outputs = FALSE)